Class 3
Zanim zaczniemy
0. Przydatne zasoby pomocowe w pracy z R - część
trzecia
- Modele liniowe w ekologii - przegląd
- Dlaczego nie musimy zagnieżdzać zmiennych?
- Równoważność metod LM oraz prostych testów statystycznych
- Modele
mieszane z
lmer() - “Statistical rethinking” - książka o modelach liniowych no. 1
- Linear models in R - wstęp
- Wszystko co należy wiedzieć o GLMM
library(tidyverse)
library(MuMIn)
library(lme4)
library(lmerTest)
library(sjPlot)
dane <- read_csv("BTdata.csv",
col_types = cols(sex = col_factor(levels = c("Fem", "Male", "UNK")),
habitat = col_factor(levels = c("forest", "park"))))
dane <- filter(dane, sex != "UNK")A. Regresja, ANOVA + ANCOVA
1. Regresja wielokrotna
We wszystkich modelach liniowych tworzonych w R
stosujemy konwencję formuły do zapisania zależności między zmiennymi w
modelu. Oto kilka podstawowych elementów formuły, wraz z ich
interpretacją:
| Formuła | Znaczenie |
|---|---|
y ~ x |
Regresja prosta (x ciągłe) lub jednokierunkowa ANOVA
(x kategoryczne); jeśli x ma 2 poziomy taka
notacja da wynik identyczny z t.test() |
y ~ x + z |
Regresja wielokrotna (x i z ciągłe);
dwukierunkowa ANOVA (x i z kategoryczne);
ANCOVA (jedna zmienna ciągła, druga kategoryczna) |
y ~ x * z |
Model liniowy z interakcją x\(\times\)y |
y ~ x + z + x:z |
Alternatywny zapis poprzedniego modelu |
y ~ (x + z + a)^2 |
Model z 3 czynnikami, i ich wszystkimi interakcjami 2 stopnia (zapis
równoważny x + z + a + x:a + x:z + z:a) |
y ~ x * z * a |
Model z 3 czynnikami i ich wszystkimi interakcjami (łącznie z
x:z:a) |
y ~ x/z |
Model ze zmienną z zagnieżdżoną w x |
y ~ x + (1|z) |
Model z efektem losowym z (notacja lmer()
oraz lme()) |
Zadanie Opierając się na prostej regresji skonstruuj
model, w którym na tarsus potencjalnie wpływają:
hatchdate, sex, habitat,
interakcja habitat\(\times\)hatchdate, interakcja
sex\(\times\)hatchdate. Zapisz go w
nowym obiekcie - czy którykolwiek czynnik modelu jest istotny
statystycznie?
Oczekiwany wynik
model1 <- lm(tarsus ~ hatchdate + habitat + sex + habitat:hatchdate + hatchdate:sex,
data = dane)
summary(model1)##
## Call:
## lm(formula = tarsus ~ hatchdate + habitat + sex + habitat:hatchdate +
## hatchdate:sex, data = dane)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.71722 -0.29330 0.05528 0.30670 2.09627
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15.231391 0.565077 26.955 <2e-16 ***
## hatchdate -0.019726 0.011697 -1.686 0.0921 .
## habitatpark -0.524629 0.649435 -0.808 0.4194
## sexMale -0.003442 0.648839 -0.005 0.9958
## hatchdate:habitatpark 0.011810 0.013430 0.879 0.3795
## hatchdate:sexMale 0.008137 0.013417 0.606 0.5444
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4631 on 775 degrees of freedom
## Multiple R-squared: 0.1561, Adjusted R-squared: 0.1506
## F-statistic: 28.66 on 5 and 775 DF, p-value: < 2.2e-16
plot(model1)Uprość model usuwając kolejno nieistotne interakcje. Jaka jest ostateczna interpretacja finalnego modelu?
Oczekiwany wynik
model1 <- lm(tarsus ~ hatchdate + habitat + sex + habitat:hatchdate,
data = dane)
summary(model1)##
## Call:
## lm(formula = tarsus ~ hatchdate + habitat + sex + habitat:hatchdate,
## data = dane)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.73620 -0.28877 0.05101 0.30354 2.09101
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15.030980 0.458192 32.805 <2e-16 ***
## hatchdate -0.015579 0.009487 -1.642 0.101
## habitatpark -0.519466 0.649114 -0.800 0.424
## sexMale 0.389543 0.033161 11.747 <2e-16 ***
## hatchdate:habitatpark 0.011705 0.013423 0.872 0.383
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4629 on 776 degrees of freedom
## Multiple R-squared: 0.1557, Adjusted R-squared: 0.1513
## F-statistic: 35.77 on 4 and 776 DF, p-value: < 2.2e-16
model1 <- lm(tarsus ~ hatchdate + habitat + sex,
data = dane)
summary(model1)##
## Call:
## lm(formula = tarsus ~ hatchdate + habitat + sex, data = dane)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.70922 -0.29145 0.05238 0.31820 2.08274
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.749240 0.324830 45.406 <2e-16 ***
## hatchdate -0.009733 0.006711 -1.450 0.147
## habitatpark 0.045817 0.033155 1.382 0.167
## sexMale 0.389645 0.033156 11.752 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4628 on 777 degrees of freedom
## Multiple R-squared: 0.1548, Adjusted R-squared: 0.1516
## F-statistic: 47.45 on 3 and 777 DF, p-value: < 2.2e-16
\(\blacksquare\)
2. Regresja wielokrotna via seleckja modeli
Aby przeprowadzić selekcję modeli musimy najpierw ponownie policzyć najbardziej złożony model - tzw. model maksymalny.
model0 <- lm(tarsus ~ hatchdate + habitat + sex + habitat:hatchdate + hatchdate:sex,
data = dane, na.action = 'na.fail')Funkcja dredge() w pakiecie MuMIn
przeprowadza automatyczną selekcję modeli, wraz z dopasowaniem modeli
zagnieżdżonych w głównym modelu stworzonym powyżej.
Zadanie Spróbuj wykonać taką selekcję, zapisując jej
wynik do obiektu model_sel. Następnie zobacz jak wygląda
powstały w ten sposób obiekt.
Oczekiwany wynik
model_sel <- dredge(model0, trace = T)## Fixed term is "(Intercept)"
## 0 : lm(formula = tarsus ~ 1, data = dane, na.action = "na.fail")
## 1 : lm(formula = tarsus ~ habitat + 1, data = dane, na.action = "na.fail")
## 2 : lm(formula = tarsus ~ hatchdate + 1, data = dane, na.action = "na.fail")
## 3 : lm(formula = tarsus ~ habitat + hatchdate + 1, data = dane, na.action = "na.fail")
## 4 : lm(formula = tarsus ~ sex + 1, data = dane, na.action = "na.fail")
## 5 : lm(formula = tarsus ~ habitat + sex + 1, data = dane, na.action = "na.fail")
## 6 : lm(formula = tarsus ~ hatchdate + sex + 1, data = dane, na.action = "na.fail")
## 7 : lm(formula = tarsus ~ habitat + hatchdate + sex + 1, data = dane,
## na.action = "na.fail")
## 11 : lm(formula = tarsus ~ habitat + hatchdate + habitat:hatchdate +
## 1, data = dane, na.action = "na.fail")
## 15 : lm(formula = tarsus ~ habitat + hatchdate + sex + habitat:hatchdate +
## 1, data = dane, na.action = "na.fail")
## 22 : lm(formula = tarsus ~ hatchdate + sex + hatchdate:sex + 1, data = dane,
## na.action = "na.fail")
## 23 : lm(formula = tarsus ~ habitat + hatchdate + sex + hatchdate:sex +
## 1, data = dane, na.action = "na.fail")
## 31 : lm(formula = tarsus ~ habitat + hatchdate + sex + habitat:hatchdate +
## hatchdate:sex + 1, data = dane, na.action = "na.fail")
model_sel## Global model call: lm(formula = tarsus ~ hatchdate + habitat + sex + habitat:hatchdate +
## hatchdate:sex, data = dane, na.action = "na.fail")
## ---
## Model selection table
## (Int) hbt htc sex hbt:htc htc:sex df logLik AICc delta weight
## 5 14.30 + 3 -506.383 1018.8 0.00 0.191
## 7 14.75 -0.009317 + 4 -505.416 1018.9 0.09 0.183
## 8 14.75 + -0.009733 + 5 -504.458 1019.0 0.20 0.173
## 6 14.28 + + 4 -505.513 1019.1 0.28 0.166
## 16 15.03 + -0.015580 + + 6 -504.075 1020.3 1.46 0.092
## 23 14.95 -0.013380 + + 5 -505.234 1020.5 1.75 0.080
## 24 14.94 + -0.013750 + + 6 -504.279 1020.7 1.87 0.075
## 32 15.23 + -0.019730 + + + 7 -503.890 1021.9 3.13 0.040
## 1 14.51 2 -570.153 1144.3 125.52 0.000
## 3 14.99 -0.010100 3 -569.187 1144.4 125.61 0.000
## 4 14.99 + -0.010530 4 -568.345 1144.7 125.94 0.000
## 2 14.48 + 3 -569.394 1144.8 126.02 0.000
## 12 15.29 + -0.016650 + 5 -567.988 1146.1 127.26 0.000
## Models ranked by AICc(x)
Jak zinterpretujesz powstałe rezultaty? Który model jest najlepszy? Czy jest to wiodące prowadzenie na tle innych modeli?
Zadanie Użyj funkcji model_avg() by
przeprowadzić uśrednianie modeli mieszczących się w zakresie deltaAICc
<= 2 (subset = delta <= 2). Czy wyniki są bardziej
jednoznaczne?
Oczekiwany wynik
model_average <- model.avg(model_sel, subset = delta <= 2)
summary(model_average)##
## Call:
## model.avg(object = model_sel, subset = delta <= 2)
##
## Component model call:
## lm(formula = tarsus ~ <7 unique rhs>, data = dane, na.action = na.fail)
##
## Component models:
## df logLik AICc delta weight
## 3 3 -506.38 1018.80 0.00 0.20
## 23 4 -505.42 1018.88 0.09 0.19
## 123 5 -504.46 1018.99 0.20 0.18
## 13 4 -505.51 1019.08 0.28 0.17
## 1234 6 -504.08 1020.26 1.46 0.10
## 235 5 -505.23 1020.55 1.75 0.08
## 1235 6 -504.28 1020.67 1.87 0.08
##
## Term codes:
## habitat hatchdate sex habitat:hatchdate hatchdate:sex
## 1 2 3 4 5
##
## Model-averaged coefficients:
## (full average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 14.638242 0.416181 0.416533 35.143 <2e-16 ***
## sexMale 0.327314 0.298596 0.298956 1.095 0.274
## hatchdate -0.007211 0.008631 0.008638 0.835 0.404
## habitatpark -0.030388 0.258184 0.258431 0.118 0.906
## habitatpark:hatchdate 0.001122 0.005398 0.005403 0.208 0.836
## hatchdate:sexMale 0.001295 0.006144 0.006151 0.211 0.833
##
## (conditional average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 14.638242 0.416181 0.416533 35.143 <2e-16 ***
## sexMale 0.327314 0.298596 0.298956 1.095 0.274
## hatchdate -0.011481 0.008341 0.008353 1.375 0.169
## habitatpark -0.057619 0.353305 0.353648 0.163 0.871
## habitatpark:hatchdate 0.011705 0.013423 0.013444 0.871 0.384
## hatchdate:sexMale 0.008035 0.013418 0.013439 0.598 0.550
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
\(\blacksquare\)
3. ANOVA
Zaimportuj do R dane pg.csv.
Zadanie Korzystając z funkcji lm()
skonstruuj model, w którym zmienną zależną jest weight a
niezależnymi group oraz mod. Jak wyglądają
wyniki takiego modelu?
Oczekiwany wynik
model3 <- lm(weight ~ group + mod, data = plantgr)
summary(model3)##
## Call:
## lm(formula = weight ~ group + mod, data = plantgr)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.9710 -0.3765 -0.0330 0.2100 1.2600
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.1410 0.2280 22.550 <2e-16 ***
## grouptrt1 -0.3710 0.2792 -1.329 0.1955
## grouptrt2 0.4940 0.2792 1.769 0.0886 .
## modN -0.2180 0.2280 -0.956 0.3478
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6244 on 26 degrees of freedom
## Multiple R-squared: 0.2891, Adjusted R-squared: 0.2071
## F-statistic: 3.525 on 3 and 26 DF, p-value: 0.02881
Zadanie Poarównaj uzyskany output z wynikiem użycia
w podobny sposób funkcji aov().
Oczekiwany wynik
model4 <- aov(weight ~ group + mod, data = plantgr)
summary(model4)## Df Sum Sq Mean Sq F value Pr(>F)
## group 2 3.766 1.8832 4.831 0.0164 *
## mod 1 0.356 0.3564 0.914 0.3478
## Residuals 26 10.136 0.3898
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Czy dodanie do modelu interakcji group:mod wyjaśnia
dodatkową zmienność modelu? Wskazówka \(R^2\) modelu policzonego za pomocą funkcji
aov() można podejrzeć stosując
summary.lm(model)$r.squared.
Oczekiwany wynik
model5 <- aov(weight ~ group * mod, data = plantgr)
summary(model5)## Df Sum Sq Mean Sq F value Pr(>F)
## group 2 3.766 1.8832 5.401 0.0116 *
## mod 1 0.356 0.3564 1.022 0.3221
## group:mod 2 1.768 0.8839 2.535 0.1003
## Residuals 24 8.368 0.3487
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary.lm(model4)$r.squared## [1] 0.2891461
summary.lm(model5)$r.squared## [1] 0.4131247
Zadanie Czy model z interakcją spełnia założenia modelu liniowego?
Oczekiwany wynik
plot(model5)\(\blacksquare\)
B. Model uogólniony
1. Regresja logistyczna - GLM
Zaimportuj dane BTEPC.txt do obiektu btepc.
Zmienne MAGEO i FAGEO powinny mieć charakter
kategoryczny (factor). Jakiego rodzaju zmienne znajdują się
w pliku?
Zadanie Skonstruuj wykres pokazujący zależność
FEPC od MTAR.
Oczekiwany wynik
plot(jitter(FEPC, 0.3) ~ MTAR, data = btepc)Zadanie Skonstruuj model testujący wpływ interakcji
długości skoku partnera socjalnego (MTAR) i wieku samicy
(FAGEO) na prawdopodobieństwo “skoku-w-bok” u sikory
modrej. Pamiętaj o zmianie rodziny rozkładów dla modelu GLM.
Spodziewany wynik
model6 <- glm(FEPC ~ MTAR*FAGEO, data = btepc, family = "binomial")
summary(model6)##
## Call:
## glm(formula = FEPC ~ MTAR * FAGEO, family = "binomial", data = btepc)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.194 -1.134 -0.595 1.197 1.925
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.9114 9.3003 -0.206 0.8372
## MTAR 0.1083 0.5448 0.199 0.8424
## FAGEO2 36.0285 17.7443 2.030 0.0423 *
## MTAR:FAGEO2 -2.1401 1.0384 -2.061 0.0393 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 149.53 on 109 degrees of freedom
## Residual deviance: 139.68 on 106 degrees of freedom
## (4 obserwacje zostały skasowane z uwagi na braki w nich zawarte)
## AIC: 147.68
##
## Number of Fisher Scoring iterations: 4
Poniższy kod pozwala zobaczyć przewidywania modelu na realnej, biologicznej skali prawdopodobieństwa.
plot_model(model6, type = "int")## Data were 'prettified'. Consider using `terms="MTAR [all]"` to get smooth plots.
\(\blacksquare\)
C. Model mieszany
1. Czy skok jest odziedziczalny?
Aby przetestować genetyczną komponentę cechy zmierzonej w ramach eksperymentu genetycznego wykorzystującego cross-fostering należy zbadań model zawierający efekty losowe gniazda wychowywania i gniazda genetycznego (rodziny, matki genetycznej). Czynniki takie typowo analizujemy jako efekty losowe - model musi więc mieć charakter mieszany.
Załaduj dane BTdata.csv, zadbaj by
fosternest oraz dam miały charakter czynników
kategorycznych (factor).
dane$fosternest <- as.factor(dane$fosternest)
dane$dam <- as.factor(dane$dam)Model dopasujemy w funkcji lmer() z pakietu
lme4 - zawierać ona będzie czynnik ustalony płci oraz
czynniki losowe fosternest i dam.
Zadanie Wiedząc, że w lmer() czynnik
losowy dodajemy wg konwencji (1|nazwa_czynnika) -
skonstruuj i policz odpowiedni model. Gdzie w jego podsumowaniu znajdują
się informacje o efektach losowych? Czy model spełnia założenia modelu
liniowego?
Spodziewany wynik
model7 <- lmer(tarsus ~ sex + (1|fosternest) + (1|dam), data = dane)
plot(model7)summary(model7)## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: tarsus ~ sex + (1 | fosternest) + (1 | dam)
## Data: dane
##
## REML criterion at convergence: 890.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.1407 -0.6128 0.0703 0.6250 4.9918
##
## Random effects:
## Groups Name Variance Std.Dev.
## dam (Intercept) 0.05401 0.2324
## fosternest (Intercept) 0.01621 0.1273
## Residual 0.14495 0.3807
## Number of obs: 781, groups: dam, 106; fosternest, 104
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 14.29877 0.03338 169.96645 428.33 <2e-16 ***
## sexMale 0.38518 0.02890 700.38893 13.33 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## sexMale -0.456
Zadanie Aby przetestować czy czynnik losowy jest
istotny statystycznie oblicz konkurencyjny model pozbawiony czynnika
rodzica genetycznego (dam) i porównaj go z wcześniejszym
modelem za pomocą funkcji anova(). Który model ma niższe
AIC i jest tym samym lepszy? Czy różnica jest istotna
statystycznie?
Spodziewany wynik
model8 <- lmer(tarsus ~ sex + (1|fosternest), data = dane)
anova(model7, model8)## refitting model(s) with ML (instead of REML)
## Data: dane
## Models:
## model8: tarsus ~ sex + (1 | fosternest)
## model7: tarsus ~ sex + (1 | fosternest) + (1 | dam)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model8 4 971.80 990.45 -481.90 963.80
## model7 5 890.44 913.75 -440.22 880.44 83.36 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
\(\blacksquare\)
2. Czy temperatura różnie wpływ na tempo wzrostu u różnych gatunków plnaktonu?
Załaduj dane PLD.txt do obiektu pld.
Przyjrzyj się strukturze danych. Zamień kolumnę species na
factor:
pld$species <- as.factor(pld$species)Zadanie Stwórz model mieszany analizujący zależność
długości okresu larwalnego (PLD) od temperatury u różnych gatunków
bezkręgowców. Użyj species jako efektu losowego,
zlogarytmuj w modelu zmienną zależną i niezależną. Czy zmienność między
gatunkami jest duża, czy mała?
Spodziewany wynik
model9 <- lmer(log(pld) ~ log(temp) + (1|species), data = pld)
summary(model9)## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: log(pld) ~ log(temp) + (1 | species)
## Data: pld
##
## REML criterion at convergence: 220.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.75675 -0.49022 -0.00731 0.47073 2.94018
##
## Random effects:
## Groups Name Variance Std.Dev.
## species (Intercept) 0.78083 0.8836
## Residual 0.04207 0.2051
## Number of obs: 214, groups: species, 72
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 6.1661 0.1746 209.0718 35.32 <2e-16 ***
## log(temp) -1.1508 0.0493 152.1852 -23.34 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## log(temp) -0.798
Zadanie Uzupełnij model możliwością zmiennych
międzygatunkowo nachyleń zależności pld~temp zamieniając
1 w definicji efektu losowego na 1+log(temp).
Który model jest lepszy?
model10 <- lmer(log(pld) ~ log(temp) + (1+log(temp)|species), data = pld)
summary(model10)## Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
## Formula: log(pld) ~ log(temp) + (1 + log(temp) | species)
## Data: pld
##
## REML criterion at convergence: 184.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.7463 -0.4046 0.0015 0.4255 3.8863
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## species (Intercept) 4.35612 2.0871
## log(temp) 0.28894 0.5375 -0.92
## Residual 0.02035 0.1426
## Number of obs: 214, groups: species, 72
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 7.07746 0.29590 58.36905 23.92 <2e-16 ***
## log(temp) -1.43319 0.08364 53.21283 -17.14 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## log(temp) -0.943
anova(model9, model10)## refitting model(s) with ML (instead of REML)
## Data: pld
## Models:
## model9: log(pld) ~ log(temp) + (1 | species)
## model10: log(pld) ~ log(temp) + (1 + log(temp) | species)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model9 4 222.02 235.48 -107.010 214.02
## model10 6 190.79 210.98 -89.394 178.79 35.232 2 2.236e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Poniższa funkcja (szczegółami zajmiemy się jutro) ilustruje badane
dane, pokazując zależności log(pld)~log(temp) dla
wszystkich obecnych w danych gatunków:
library(ggplot2)
qplot(x = log(temp), y = log(pld),
data = pld, col = species, ) +
geom_smooth(method = "lm", se = F) +
theme_bw() + theme(legend.position="none")## `geom_smooth()` using formula 'y ~ x'
\(\blacksquare\)
D. Expressem przez bootstrapping
Do tej pory do naszych analiz wykorzystywaliśmy metody parametryczne
- możemy jednak spróbować wykorzystać inną, prostszą metodę, która
pozwoli nam uwolnić się od ograniczeń rozkładowych naszych danych.
Spójrzmy na przykład na rozkład zmiennej hatchdate - nie
jest on najpiękniejszy i na pierwszy rzut oka nie przypomina on
zgrabnego rozkładu normalnego:
hist(dane$hatchdate, 15)Spróbujmy skonstruować test oparty na bootstrappingu, którym
sprawdzimy czy ptaki z lasu (forest) różnią się pod kątem
daty klucia od ptaków z parków (park). Szybki rzut oka na
wykres pudełkowy nie wskazuje na zbyt wielkie różnice.
plot(hatchdate ~ habitat, data = dane)- Do skonstruowania bootstrappingu wykorzystamy pętlę
for, która pozwala powtarzać jakąś operację wiele razy. Jeśli chcielibyśmy np. wydrukować na ekranie (print()) jakiś tekst pięć razy - moglibyśmy osiągnąć to w następujący sposób:
for(i in 1:5) {
print("Tekst najbardziej oryginalny!")
}## [1] "Tekst najbardziej oryginalny!"
## [1] "Tekst najbardziej oryginalny!"
## [1] "Tekst najbardziej oryginalny!"
## [1] "Tekst najbardziej oryginalny!"
## [1] "Tekst najbardziej oryginalny!"
Zauważ, że pętla for wykorzystuje jakąś zmienną (tutaj
i), zmieniając jej wartość w zadanym zakresie (tutaj
1:5) i za każdym razem (dla każ∂ego i)
wykonując jakieś działania ujęte wewnątrz {}.
- Potrzebna nam będzie również metoda losowego wybierania pewnych
wierszy z naszych oryginalnych danych - ze zwracaniem. Idealnie do tego
nadaje się funkcja
sample(). Zobaczmy, jak losuje ona zadaną liczbę wartości (10) z podanej listy (liczby 1 do 20), ze zwracaniem (replace = T):
sample(1:20, size = 10, replace = T)## [1] 17 20 12 16 2 17 13 17 18 18
Uwaga! Funkcja ta używa generatora liczb pseudolosowych - więc u każdego wynik tej operacji będzie inny!
- Wreszcie - musimy mieć naszą statystykę testową. Niech będzie nią różnica między średnią datą klucia w sikor leśnych i parkowych \(d_0\) - jeśli wynosi ona zero, ptaki mają identyczne daty klucia:
d0 <- mean(filter(dane, habitat == "forest")$hatchdate) -
mean(filter(dane, habitat == "park")$hatchdate)
d0## [1] -0.2214331
Ujemna różnica wskazuje na nieco wyższą wartość (późniejsze klucie) w parku - pytanie, czy istotnie późniejsze?
Oto co powinna robić nasza procedura bootstrappingu:
- Krok 1: wylosuj z zestawu danych k = 781 obserwacji (wierszy) ze zwracaniem (tyle wierszy składa się na oryginalny zestaw danych).
- Krok 2: dla każdego z tych wylosowanych zestawów policz jego własną
statystykę
d- różnicę w średniej dacie klucia między lasem i parkiem. - Krok 3: Powtórz całą procedurę N = 1000 razy (jest to arbitralna
wielkość, im więcej - tym lepiej). Za każdym razem zapisz wynik do
zbiorczego obiekty
out. - Krok 4: Skonstruuj rozkład statystyki
di zobacz gdzie wypada w niej jej oryginalna (oparta na wejściowych danych) wartośćd0.
N <- 1000 # liczba zakładanych powtórzeń
k <- nrow(dane) # wielkość oryginalnego zbioru danych
# oryginalna statystyka testowa
d0 <- mean(filter(dane, habitat == "forest")$hatchdate) -
mean(filter(dane, habitat == "park")$hatchdate)
# obiekt wynikowy przygotowany na N=1000 końcowych wartości
out <- numeric(N)
for(i in 1:N) {
# losujemy 828 wierszy, ze zwracaniem
# a oryginalnej puli unikatowych 828 wierszy danych
rows <- sample(1:k, size = k, replace = T)
# tworzymy tymczasowe dane wewnątrz pętli
# zawierające tylko wylosowane wiersze
temp_dane <- dane[rows, ]
# oblicza na podstawie tak przygotowanych danych
# tymczasową statystykę testową dla powtórzenia
# nr i pętli
d_temp <- mean(filter(temp_dane, habitat == "forest")$hatchdate) -
mean(filter(temp_dane, habitat == "park")$hatchdate)
# zapisujemy wyliczoną statystykę
# do zmiennej out w miejscu i
out[i] <- d_temp
}Wyświetlmy rozkład symulowanej statystyki d oraz
położenie d0 na rozkładzie:
hist(out, 50)
abline(v = d0, lwd = 2)
abline(v = 0, col = "red", lwd = 2, lty = 2)Jak podejdziesz to odpowiedzi na pytanie - czy ta symulacja wskazuje na istotną statystycznie różnicę między symulowaną statystyką a wartością zero?
\(\blacksquare\)